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We consider the grand potential f2 of a two-dimensional weakly interacting homogeneous Bose 
gas at zero temperature. Building on a number-conserving Bogoliubov method for a lattice model 
in the grand canonical ensemble, we calculate the next order term as compared to the Bogoliubov 
prediction, in a systematic expansion of Q in powers of the parameter measuring the weakness of 
the interaction. Our prediction is in very good agreement with recent Monte Carlo calculations. 
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Recent experimental progress with ultracold atoms has 
renewed the interest in the two-dimensional weakly inter- 
acting Bose gas [TJ [5] . In view of a comparison to future 
experimental results on the equation of state of the gas at 
low temperatures, this raises the question of the accuracy 
of existing theoretical work [3j |4] . Since the pioneering 
works of Schick [S] and Popov [5] on the energy of the 
weakly interacting Bose gas in two dimensions, several 
recent predictions have been obtained. In mathemati- 
cal physics, it was proved that Schick's formula for the 
ground state energy is asymptotically exact in the limit 
of vanishing density [7] . Numerically, very precise Monte 
Carlo calculations of the ground state energy have been 
performed [SlIH]. Analytically, Popov's result was con- 
firmed by a Bogoliubov type theory [TU] (and by Monte 
Carlo calculations |9]) but several attempts to calculate 
analytically the energy beyond Popov's result have led to 
non-identical predictions [TTJ [TU |T3j . 

The most systematic among the theoretical approaches 
are those relying on an expansion of the energy in powers 
of a small parameter. This is the case of the approaches 
[HIlO], which have led to the equation of state [14] 
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where m is the mass of a boson, a > is the two- 
dimensional scattering length among the particles, fi is 
the chemical potential and 7 = 0.57721566 ... is Eu- 
ler's constant. Remarkably Eq.([lJ is universal, i.e. it 
depends on the interaction potential through the scat- 
tering length only. One obtains from ([ij the grand po- 
tential ^ = E — fiN in the thermodynamic limit, where 
E is the gas energy and N the atom number, by a simple 
integration over /i since N 
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where is the surface of the gas. The small parameter 
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is apparent in ([2|. 

In the vanishing density (or chemical potential) limit, 
the prediction ([T]) can be checked to be asymptotically 
equivalent to Schick's formula, as it should be. In a fur- 
ther expansion of the energy in terms of the density, the 
precise value of the constant under the logarithm in ([T| 
eventually matters. In particular, a careful account of 
the low-energy two body T-matrix is essential to derive 
this constant [T3]. Its value turns out to agree with re- 
cent Monte Carlo results [3], it however differs from the 
prediction of Ref . [12] . For not extremely small values of 
the density, a significant deviation is observed between 
the energy deduced from Eq. ([2| and Monte Carlo re- 
sults [HI in); see in Fig|l]the fact that the symbols signif- 
icantly deviate from unity. Furthermore, this deviation 
is not accounted for by the beyond-Bogoliubov theories 
of Refs. [m [13], see in Figjljthe fact that the symbols 
significantly deviate from the dashed and dotted lines. 

In the present work we extend the Bogoliubov 
method [TO] as in [T5] and we go one step further than 
Eq. ([2| in the expansion of the grand potential in powers 
of e(/i). We obtain in the thermodynamic limit: 
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where the numerical constant / is given by a multiple 
integral that we have evaluated numerically: 
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Since the extra term that we have found with respect to 
([2]) is indeed o(l), this analytically confirms that the nu- 
merical constant inside the logarithm in ([iJ is the correct 
one. Furthermore, the inclusion of the extra term leads 
to a now satisfactory agreement with the numerical re- 
sults of [Hl[n]) see in Figjljthe agreement of the solid line 
with the plotting symbols. 

Our model: In a first stage we still consider a general 
value d of the space dimension. As a regularization 
scheme to treat ultraviolet divergences, we use a lattice 
model [TU] to represent the interacting Bose gas, with the 
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FIG. 1; Ground state energy i5 of a two dimensional Bose gas, 
as a function of the gas density, and in units of the Bogoli- 
ubov prediction iSsog resulting from Solid line: energy 
obtained from the beyond Bogoliubov analytical prediction 
Q derived in this work. Dashed line: analytical predic- 
tion of [TT] taking the exact expression f{u) = u + + 
2u'^e^^^Fii{~2/u) of the function / introduced in QT], where 
Ei is the function exponential integral. Dotted line: ana- 
lytical prediction of [13] in the form [16]. Plotting symbols 
with error bars: numerical results of [8], for interactions given 
by hard disks (crosses) and by soft disks (circles); numerical 
results of [5], for dipolar interactions (diamonds). We have 
restricted the results of [8] to their universal range, that is to 
low enough values of pa^ such that the hard disks and soft 
disks models give the same values of the energy within the er- 
ror bars. For the same reason, we have restricted the results 
of to pa^ < 10~®. The inset is a magnification. 
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where the limit £ — > is eventually taken to recover the 
continuous gas. The discrete positions r run over the 
lattice ^Z'^, where £ is the lattice spacing; there is no 
trapping potential, but a quantization box [0,L]'^ with 
periodic boundary conditions. The bosonic field obeys 
the discrete commutation relations [^(ri), ■i/;'''(r2)] — 
'^n.ra/^''- In the kinetic energy operator, A is a discrete 
representation of the Laplacian on the lattice, such that 
the plane wave e''^ '' is an eigenstate of A with the eigen- 
value —k'^, the discrete nature of the lattice allowing one 
to restrict the values of k to the first Brillouin zone 
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The on-site interactions are characterized by a coupling 
constant go adjusted to reproduce the correct value of the 
scattering length a in the two-body scattering problem, 
as detailed in [T7]. For the two-dimensional case one 
obtains 
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We shall find in 2D that the first correction beyond ^ is 
universal, it depends on the interaction only through the 
scattering length a in the zero chemical potential limit. 
The microscopic details of our model, the fact that it is a 
lattice model or that Catalan's constant G — 0.91596 . . . 
appears in ([s]), are thus not relevant in this limit. 
Elimination of the condensate mode: We now assume 
that the ground state of H in the thermodynamic limit 
is a condensate, so that we take d=2[51[5]orc? = 3. 
We then use Bogoliubov method to eliminate the con- 
densate mode and obtain a Hamiltonian for the field of 
non-condensed particles. We use here a U{1) symmetry 
preserving approach in the spirit of [TSt [19] , adjusted to 
the case of the grand canonical ensemble. We split the 
field operator as the sum of the condensate field and the 
field of the non-condensed modes: 



V>(r) = 0(r)ao + V'nc(r), 



(9) 



where ao is the annihilation operator in the condensate 
mode cj){r) = 1/L''/^. We eliminate the condensate par- 
ticle number nk^o = o-qCio using 



"-k=o ^ N - iVnc, 
where N is the total number of particles and 
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is the number of non-condensed particles. To complete 
the elimination of the condensate mode, we introduce the 
representation [20] 



0.0 = A nJ^/Jo with A={1 + nk=o)~^/^ 



ao. 
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As shown in Eq.(5.40) of [20 one has the exact relations 
AA^ = 1 and A''' A = 1 — |vaco)(vaco|, where |vaco) is the 
vacuum state for the condensate mode. The condensate 
mode elimination is completed by inclusion of A^ in the 
non-condensed field, defining as in |l9j the field operator 
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which conserves the total particle number. In the ther- 
modynamic limit, the condensate mode has a vanishing 
probability to be empty, so that A^A AA^ = 1, and A 
and A^ obey simple commutation relations in this limit. 

In the canonical ensemble, it remains to inject the split- 
ting ([9| into H and to eliminate the condensate mode, 
finally replacing the operator N by its known value N. 
One obtains contributions of various degrees in A, start- 
ing from degree two. The terms of degree two in A gives 
the Bogoliubov Hamiltonian, the terms of higher degrees 
may be treated by perturbation theory. 

In the grand canonical ensemble, however, the chemical 
potential fi rather than the particle number is known; 
at zero temperature, N does not fluctuate but assumes 
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an a priori unknown value N{fj,), a function of n to be 
determined order by order in the weakly interacting limit. 
To zeroth order in A, the gas is a pure condensate and 
one obtains the mean-field type relation 
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It is then convenient to split N as 

N{^i) = M'^\^I) + ^N{^I). 
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As we shall see, (5-/V(/i) to leading order is second order in 
A, as the mean number of non-condensed particles (iVnc). 

After some calculation, neglecting unity as compared 
to the condensate atom number and replacing N with 
N(fi), we obtain the desired rewriting of the Hamiltonian 
with no reference to the condensate mode: 
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We have used the fact that, in the spatially homogeneous 
case, one exactly has X]r^''A(r) — 0, and we have set 
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Perturbative expansion: We now expand (16 1 in powers 



of A. Keeping terms up to second order in A we obtain 



i7<2 = -i/.iV(0)(M)-H^. 
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This plays the role of the Bogoliubov Hamiltonian in the 
usual theory. Its ground state energy we thus call the 
Bogoliubov approximation for the grand potential: 
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where we have replaced N^^'>{fj,) by its value and we have 
introduced the Bogoliubov modal amplitudes obeying 
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and the corresponding Bogoliubov energies 
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Taking minus the derivative of (19) with respect to fj, to 



(14) obtain the atom number, and using 



^^.iekV^) = -VkiUk + Vk), 
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one recovers, in the thermodynamic limit, Eq.(152) in 
[TD], and thus (0 in the limit £^0 

To go beyond Bogoliubov, we collect into H3 the terms 
of degree three in A and into H4 the terms of degree four 
in A, keeping in mind that (5-/V(/i) is to leading order of 
degree two, SN{n) = TV^^)^^) _|_ . . go that 



90 

9o 
2L^ 

It 



A^(o)(Ai) 



1/2 



7V(2)(/.) 



^£''At(A + At)A 



iv(2)(^)-4 



| + |^£rfAtAtAA. 



(23) 
(24) 



We treat H4 to first order in perturbation theory and 
H3 to second order, to obtain the first correction to the 
Bogoliubov prediction for the grand potential: 

Sn - {H,) + {H,- ^-^H,) (25) 

"Bog — -n<2 

where the expectation value is taken in the ground state 
of H<2, that is in the vacuum of the operators &k appear- 
ing in the modal expansion 



A(r) = L-^^^ ^kf/fce''^ '- -I- fot y^g- 
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To find the value of TV^^) (^) ^ 

we minimize {H4) over 
keeping in mind that i?<2 and do not depend on N^'^\ 
In the thermodynamic limit, one has (N^c) — (Mic)^ and 

(7Vnel)ci(iV„e)(X),SO that 
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= - Vk{Uk + Vk). 
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The divergence of N'^^^fj,) when £ ^ is removed in the 
combination N^^\fj,) + A^(^)(^). The resulting density is 
in agreement with Eq. ([ij . One is then left with 
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Using the modal expansion (|26| and Wick's theorem we 
finally obtain after some calculation 
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We have introduced the vector k2 G 2? such that ki + 
k2 + k3 e {2tt/1)Z'^. The notation Ui, i G {1,2,3}, stands 
for Uki ■ The sum over a runs over the permutation group 
S-i. For convenience we have added iV(2)(^) to N'-°\h) 
in the denominator of the overall factor in ( p9| ) , which is 
allowed at the present order of the calculation. 
Absence of divergences in 2D: We now take the thermo- 
dynamic limit, replacing sums over k by integrals over 



the domain 2? in ( 29 1 . We also take the zero lattice spac- 
ing limit £ [22\ so that the integration domain over 
k is now IR''. Since k2 = — (ki -t- ks) and the integrand 
depends only on the moduli ki,k2 and fcs, see (29), we 



are left with a triple integral over fci , and the angle be- 
tween the vectors ki and ka. In 2D, we show below that 
this integral converges, that is it has neither an infrared 
nor an ultraviolet divergence. The first correction beyond 
the Bogoliubov energy is thus universal in 2D. Since con- 
vergence is established, we can resort to numerical inte- 
gration. After the change of variables qi — hki/{2miiY^'^ 
and pulling out a factor n[{2mfiy/^ L/{2nh)]'^ / fi, we get 
([5|. Summing r^Bog to we then obtain (|4|. 

To show the infrared convergence, we replace the inte- 
grand by its leading low-fc^ behavior: Ui and Vi diverge as 
l/\/~ki, Si vanishes as ^/kl and vanishes as ki. Including 
the Jacobian factors ki and k^ from 2D integration in po- 
lar coordinates, we see that the integral of the first term 



in the curly brackets of (29 1 converges. The contribution 



of the term due to permutation a scales as 
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so its integral over ki and k^ is also convergent. 

For the ultraviolet convergence, the full reasoning is 
rather long so we give a simplified explanation. We 



approximate each term in the integrand in ( 29 1 by its 
leading high-fc^ behavior, Ui and Si tending to unity, Vi 
vanishing as 1/kf and diverging as fcf. In the sum 
over cr, the terms with (7(3) ^ 3 are not dangerous. E.g. 
for (t(3) = 1, a factor ViVij appears, and including the 
Jacobian factors, one obtains a contribution scaling as 

'^'^ • ' <:^, (31) 



kl 



so that the resulting double integral over fci and k^ is 
convergent at infinity. The dangerous terms in the sum 
over a thus correspond to cr(3) = 3: the factor ensures 
convergence of the integral over ka over a ki-independent 
range ^ {mn/h^Y^'^. At large ki, the energy denom- 
inator ci + £2 + £3 approaches 2ei. Then, from the 
asymptotic relation Vi ~ — /i/(2ei), we see that the two 
dangerous contributions coming from the permutations 
with cr(3) — 3 exactly compensate with the first term 
~ 2V^Vi/n in the curly brackets of (29 1, which avoids 
an ultraviolet divergence of 50.. 

In conclusion, we have calculated analytically and in 
a systematic way the first correction to the Bogoliubov 



prediction for the ground state grand potential of a 2D 
weakly interacting Bose gas. We find that this correc- 
tion is universal, depending on the interaction potential 
through the scattering length only. It allows to describe 
analytically the not extremely weakly interacting regime, 
and contrarily to other analytical works, we obtain a pre- 
diction for the ground state energy in excellent agreement 
with the numerical results of [H [9] over the range where 
the results of [HI [5] are model independent. 

This work was stimulated by discussions with Elliot 
Lieb and Jakob Yngvason. The group of Y.C. is a mem- 
ber of IFRAF. 



[1] Z. Hadzibabic et al, Nature 441, 1118 (2006). 

[2] P. Clade et al, arXiv:0805.3519 (2008). 

[3] A. Posazhennikova, Rev. Mod. Phys. 78, 1111 (2006). 

[4] I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 

80, 885 (2008). 
[5] M. Schick, Phys. Rev. A 3, 1067 (1971). 
[6] V.N. Popov, Theor. Math. Phys. 11, 565 (1972). 
[7] E. Lieb and J. Yngvason, J. Stat. Phys. 103, 509 (2001). 
[8] S. Pilati, J. Boronat, J. Casulleras and S. Giorgini, Phys. 

Rev. A 71, 023605 (2005). 
[9] G.E. Astrakharchik et al, arXiv:0812.3844i^l (2008). 
[10] C. Mora and Y. Castin, Phys. Rev. A 67, 053615 (2003). 
[11] A.Yu. Cherny and A.A. Shanenko, Phys. Rev. E 64, 

027105 (2001). 
[12] J.O. Andersen, Eur. Phys. J. B 28, 389 (2002). 
[13] L. Pricoupenko, Phys. Rev. A 70, 013601 (2004). 
[14] We found that the quantity eo appearing in Popov's the- 
ory is simply eo = 4fi^/(mo^e^'^). This is obtained by 
comparing the low-fc expression of the two-body T ma- 
trix in Eq.(3.2) of [6^ with Eq.(160) of [TD]. 
[15] T. T. Wu, Phys. Rev. 115, 1390 (1959). 
[16] We used the relations nifia jh^ =t\\- l/ln(gt)] and 
= t^[-ln((jt) 4- 1/2] with q = exp(7)/2 and < 
t <l/q (L. Pricoupenko, private communication). 
[17] Y. Castin, J. Phys. IV (France) 116, 89 (2004). 
[18] C.W. Gardiner, Phys. Rev. A 56, 1414 (1997). 
[19] Y. Castin, R. Dum, Phys. Rev. A 57, 3008 (1998). 
[20] P. Carruthers, M. Nieto, Rev. Mod. Phys. 40, 411 (1968). 
[21] Setting qi = hki/(2mfi)^^^ , we used e.g. I2S2V2I = |si — 
1| < 2/(2 + gi) and (ei - e2)/(Q? - d) < 1 -h 2/(gf g|). 
[22] More precisely, in 2D, our expansion in powers of A 
is an expansion in powers of e(/i) for a fixed value of 
T]{fi) = ln{^/e) > 1, where h^/{m^'^) = /i. The Bogoli- 
ubov method indeed rehes on a Born expansion [23j (here 
in powers of go) of the scattering amplitude fk for k ~ 
1/5. From Eq.(167) of [10] the small parameter for the 
Born expansion is 2rie. This small parameter is explicitly 
obtained in our approac h, fr om the requirement A'''-^-' <^ 
iV(°): Estimating from ^ Y<^' ~ riL'^mfi/{2nh'^) for 
^ < ^, we get iV(2)/jv{0) ^ 2r?e/(l + 2r;e). One then takes 
the limit -q +00 in each coefficient of the expansion 
of Q in powers of e. This limit is exponentially fast ap- 
proached in 77: we find that = 7 is more than large 
enough, it gives the value of J in ([H]) at the 10~^ level. 
[23] E.H. Lieb and W. Liniger, Phys. Rev. 130, 1605 (1963). 



